/* Creates the data to reproduce Figure 1. 
   Run only after MONTE.PRG has been run.
   Jason Wittenberg, 12/20/2001, witty@polisci.wisc.edu
*/

new;

loadm ds=fullmc;                                @ Load data matrix created in MONTE.PRG @
v1=ds[.,1];                                     @ Assign variables @
v2=ds[.,2];
v3=ds[.,3];
v1hatkk=ds[.,4:1003];
v2hatkk=ds[.,1004:2003];
v3hatkk=ds[.,2004:3003];
v1hatsur=ds[.,3004:4003];
v2hatsur=ds[.,4004:5003];
v3hatsur=ds[.,5004:6003];

v1s={};                                         @ Initialize useful variables @
v2s={};
v3s={};
v1k={};
v2k={};
v3k={};
indx={};


trange=0.25;
do while trange <= 5;                           @ Loop over ranges @
  "Range= " trange;
    { t1s } = vtebias(v1hatsur,trange,v1);      @ How close is SUR for party 1? @
    { t2s } = vtebias(v2hatsur,trange,v2);
    { t3s } = vtebias(v3hatsur,trange,v3);
    { t1k } = vtebias(v1hatkk,trange,v1);       @ How close is KK for party 1? @
    { t2k } = vtebias(v2hatkk,trange,v2);
    { t3k } = vtebias(v3hatkk,trange,v3);
    indx=indx | trange;
    v1s=v1s~t1s;                                @ Pool the results @
    v2s=v2s~t2s;    
    v3s=v3s~t3s;
    v1k=v1k~t1k;
    v2k=v2k~t2k;
    v3k=v3k~t3k;
    trange=trange+0.25;
endo;

v1s=meanc(v1s);                                 @ Average across districts @
v2s=meanc(v2s);
v3s=meanc(v3s);

v1k=meanc(v1k);
v2k=meanc(v2k);
v3k=meanc(v3k);

screen off;
output file=figure1.out reset;                  @ Create output file for graphing @
  indx~v1s~v2s~v3s~v1k~v2k~v3k;
output off;
screen on;

"Data written to figure1.out";


/************************************************************
Jason Wittenberg, 12/20/2001, witty@polisci.wisc.edu
Inputs: 
        POST, a matrix of posterior values (N x MSIMS)
        RANGE, range (in %) around truth
        TRUTH, vector of true values
Outputs:
        VFRAC - Frac of posterior within +- range of truth
*************************************************************/
proc (1) = vtebias(post,range,truth);
   local lb,ub,vfrac;
   lb=truth-range;
   ub=truth+range;
   vfrac=(post .>= lb .and post .<= ub);
   vfrac=vfrac';
   vfrac=sumc(vfrac)./rows(vfrac);  
   retp(vfrac);                             
endp;
